Kinetics and mechanism of proton transport across membrane nanopores 
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We use computer simulations to study the kinetics and mechanism of proton passage through 
a narrow-pore carbon-nanotube membrane separating reservoirs of liquid water. Free energy and 
rate constant calculations show that protons move across the membrane diffusively in single-file 
chains of hydrogen-bonded water molecules. Proton passage through the membrane is opposed by 
a high barrier along the effective potential, reflecting the large electrostatic penalty for desolvation 
and reminiscent of charge exclusion in biological water channels. At neutral pH, we estimate a 
translocation rate of about 1 proton per hour and tube. 



Long-range proton transfer is central to processes as di- 
verse as hydrogen fuel cells [l|, Q , the enzymatic function 
of many proteins, and in particular membrane biophysics 
To explore the fundamental question of water- 
mediated proton transfer, and to design robust proton 
conducting media for technological applications, study- 
ing simpler model systems is essential. The quasi-one- 
dimensional water chains forming inside carbon nano- 
tubes Q have attracted considerable attention, with 
computer simulations suggestingproton mobilities ex- 
ceeding those even of bulk water [s, @, 0, 0, 0, ■ How- 
ever, large conductivity requires in addition a high den- 
sity of charge carriers, which depends on the free energy 
penalty required to remove protons from the bulk liquid 
and introduce them into the pores. This then raises the 
question if water-filled nanotubes can actually carry pro- 
tonic currents of high density, i.e., whether the electro- 
static desolvation penalty of the proton is compensated, 
at least in part, by its exceptionally high mobility. 

Here, we will use computer simulations to explore the 
kinetics and mechanism of proton translocation through 
nanopores. In our simulations, four rigid (6,6) armchair- 
type carbon nanotubes of 144 carbon atoms each are 
packed into a hexagonal array to form a nanotube mem- 
brane in the periodically replicated simulation box (Fig. 
[1]) . The size of the simulation box in the ^-direction par- 
allel to the tube axes is 34.3 A, and 22.5 A and 19.5 A, 
respectively, in the x- and y-directions. The membrane 
is immersed in a bath of 292 water molecules containing 
one excess proton. At T = 300 K and a density cor- 
responding to that of liquid water, the ~8-A diameter 
pores fill with single-file chains of six hydrogen bonded 
water molecules. In our simulations, the equations of 
motion are integrated with the velocity Verlet algorithm 
using a time step of 0.489 fs and a hydrogen mass of 2 




FIG. 1: (color online) . Top: Side view of the carbon nanotube 
membrane immersed in liquid water. One carbon nanotube 
is cut open to expose the chain of hydrogen bonded water 
molecules traversing the pore. Bottom: Enlarged view of a 
typical configuration of a protonic defect in the water chain 
inside the pore. 



a.m.u. For the interactions of the water molecules and 
the excess proton, we use the multistate empirical valence 
bond model (EVB) developed by Voth and collaborators 
ll| based on prior work of Warshel JJj. This model 
accurately describes the energetics of bond breaking and 
formation during aqueous proton transfer and is compu- 
tationally far less expensive than ah initio methodologies 
[sj. The water oxygen atoms interact with the carbon 
atoms of the nanotube through a Lennard-Jones poten- 
tial with e = 0.1143 kcal/mol and a — 3.27 A yielding a 
channel aperture of approximately 2 A. Periodic bound- 
ary conditions with Ewald sums for the Coulombic in- 
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teractions apply in all three spatial directions. We stress 
that the model and setup used here does not bias the sim- 
ulation towards a particular H+ transport mechanism. 

In the bulk liquid outside the carbon nanotube mem- 
brane the excess proton moves primarily as a high mobil- 
ity charge defect by proton transfers along the hydrogen 
bond network percolating through the liquid. During 
this so-called Grotthuss-process, the hydrated proton ex- 
ists in a continuum of structures including as the limiting 
cases the Eigen cation H9O4', consisting of a hydronium 
ion H3O"'" tightly hydrogen-bonded to three neighboring 
water molecules, and the Zundel cation H5OJ, in which 
the excess proton is shared between two water molecules 
[13L [3, This structural diffusion process is rapid 

with typical proton hopping times on the picosecond 
timescale such that during the nanosecond simulations 
of this study each water molecule outside the membrane 
is visited several times by the excess charge. During these 
simulations, however, the proton never entered the mem- 
brane pores. 

To clarify what prevents the proton from penetrating 
into the membrane interior despite the high proton mo- 
bility along one-dimensional water chains [8| we have cal- 
culated the free energy profile F{z) for the excess charge 
as a function of its position z along the tube axis as shown 
in Fig. [2l We computed the free energy F{z) inside the 
pore (|z| < 7.4 A) using umbrella sampling Monte Carlo 
simulations in 10 separate windows. New co nfig urations 
were generated with path sampling moves 16, |T7l|. In 
each window a total of 30,000 path shooting and shifting 
moves of 14.6 fs long trajectory segments were carried 
out amounting to a total simulation time of about 350 ps 
per window. Within the windows the free energy F{z) 
was determined from the histogram P{z) of the position 
of the center of charge [l^ , essentially the position of the 
hydronium ion averaged over all EVB states. The overall 
free energy profile was obtained by matching the free en- 
ergies calculated in the separate windows and the 2.2 ns 
equilibrium run. 

Coming from the bulk, the free energy F{z) first de- 
creases, goes through a minimum at \z\ « 10 A, and then 
rises almost linearly, reaching an approximately flat and 
8 A-wide top near the tube center. The total free en- 
ergetic cost of moving the excess charge from the bulk 
phase to the tube center is ^10 kcal/mol, about 1/3 the 
cost for a sodium ion in a similar system [l^. As the 
motion of a proton along an isolated hydrogen bonded 
water chain is an essentially barrier-less process, this free 
energy penalty is due to the desolvation energy required 
to extract the proton from the favorable bulk environ- 
ment and move it into the less polar interior of the pore, 
where the excess charge is coordinated by only two wa- 
ter molecules. Note, however, that the effective charge 
of the proton and hence also its desolvation penalty are 
substantially reduced by the dipolar polarization of the 
water chain Isl, as discussed below. This effect is absent 
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FIG. 2: Top: (color online). Free energy profile for the cen- 
ter of charge along z, where z = corresponds to the tube 
center. Bottom: Corresponding probability distribution of 
the excess charge (solid line) and a particular water oxygen 
(dashed line). The inset shows the a;j/-projected probabil- 
ity density of the protonic defect in a slab 7A < |2:| < 11 A 
just outside the carbon nanotube membrane, obtained from 
a molecular dynamics simulation of 2 ns. Red color indicates 
low proton density and blue high proton density. The rims of 
the carbon nanotubes are visible as red circles. 



for other ionic species such as the sodium ion of Ref. |l9|. 

Whereas continuum electrostatics predicts the most fa- 
vorable position of the excess charge to be deep within 
the bulk liquid, the proton has an enhanced probability 
to be located near the apolar membrane (Fig. [5] bot- 
tom). That the solvated proton is preferentially located 
near interfaces has been observed earlier in simulations 
[20I and is consistent with experiments 21 1. As depicted 
in the inset of Fig. [21 the excess charge appears to be 
located predominantly in two positions: either at the en- 
trance of the C-nanotube or in the spaces between the 
nanotubes. At both positions, the proton exists in its 
preferred Eigen-like configuration, in which the central 
hydronium ion donates three hydrogen bonds to water 
molecules, but accepts none. Such configurations occur 
also in the bulk liquid 15| , but there they are less stable 
as they strain the hydrogen bond pattern of the surround- 
ing liquid. 

To study the mechanism and kinetics of the proton 
translocation process in detail we have carried out a rate 
constant calculation using the reactive flux approach of 
Bennett and Chandler 



22l . |23[ . As a reaction coordinate 



we chose the position of the center of charge along the 
tube axis and we placed the dividing surface at the tube 
center perpendicular to its axis. A total of 5000 trajec- 
tories were initiated from initial conditions generated in 
a molecular dynamics simulation with a parabolic bias 
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FIG. 3: Reactive flux k{t) and transmission coefficient (inset). 



that kept the z-coordinate of the center of charge near 
the barrier top. The forces resulting from the bias on the 
center of charge were calculated with first order pertur- 
An uncorrelated subset of the con- 
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bation theory 

figuration with z = was then used as initial conditions 
for the reactive flux calculation. Initial momenta were 
drawn from an appropriate Maxwell-Boltzmann distri- 
bution. Trajectories were terminated at |z| = 8 A, from 
where the probability of return to the dividing surface is 
negligible. 

The reactive flux k{t) calculated from 5000 trajectories 
is shown in Fig. [31 The plateau value of k{t) is the trans- 
mission rate constant k w 6.4 x lO^s^^ for one tube. For 
a proton concentration corresponding to pH=7 one ob- 
tains a protonic current of about 1 proton per hour and 
tube. We find that proton transport through the nano- 
tube membrane is positively correlated with water flow. 
The corresponding electro-osmotic drag coefficient ifdrag 
[l[ is between about 0.5 and 1 water molecule per trans- 
ported proton, as estimated from the correlated displace- 
ments of the proton and water chain in the reactive-flux 
simulations. 

From the calculated rate k of directional proton 
translocations per pore in the absence of electric fields, 
one can estimate the proton conductivity a of (6,6) nano- 
tube membranes. In the linear response limit, the num- 
ber of protons translocated per pore is « kfieV where 
V is the applied voltage and (3 = l/kgT. For an area 
density p « 10^* m"^ of nanotubes, the current density 
becomes pkfSe^ « 100 A m^^V~^ at room temperature 
for a rate fc « 15 s~^ for pH«2. This is about two orders 
of magnitude below those of polymer electrolyte mem- 
branes used in fuel cells However, this estimate ig- 
nores that the rate of proton translocation should here 
grow exponentially with applied voltage, as it is deter- 
mined largely by proton desolvation (i.e., the low charge 
carrier concentration in the membrane) and not by the 
high proton mobility in the nanotubes. 

As the protonic defect passes through the pore, it ef- 
fectively fiips the dipolar orientation of the water chain. 



This dipole inversion is associated with a displacement 
current traveling in the direction opposite to the proton 
motion. As a consequence, the efii'ective charge trans- 
ported through the membrane by proton translocation 
alone is only about ~60% of an elementary charge [1]. 
Proton transfer across the membrane is completed when 
the orientation of the original dipole chain is restored 
by a hydrogen bonding defect passing through the pore 
[5|. This hydrogen bonding defect carries the remaining 
^^40% of the elementary charge and its passage prepares 
the water chain for transport of the next proton. In sep- 
arate simulations of a system of 4 nanotubes and 292 
TIP3P water molecules, we observed three reorientations 
during 15 ns, corresponding to a rate of about 1/(20 ns) 
per tube. This is considerably slower than the rate of 
dipolar reorientation in isolated tubes, '--^1/(2 ns) [1,[2^, 
refiecting the fact that reorientation proceeds through 
movement of a hydrogen-bond defect that carries an ef- 
fective charge through the low-polarity membrane. How- 
ever, dipole reorientation is still much faster than proton 
transfer, and thus not rate limiting. 

The rather low transmission coefficient of k « 0.065 
(inset of Fig. [3]) found in our simulations may origi- 
nate from two different causes. Either the position of 
the center of charge is not a suitable reaction coordinate 
capable of capturing the essential transition mechanism 
or the transition is of diffusive nature [l3| . In both cases 
frequent recrossings of the dividing surface reduce the 
transmission coefficient, albeit for very different reasons. 
We can distinguish these two cases by analyzing the tra- 
jectories started from the dividing surface. These tra- 
jectories were generated in pairs starting from the same 
configuration with momenta of identical magnitudes but 
opposite directions. In 54% of all pairs the two trajecto- 
ries reached different sides of the membrane and in 46% 
both trajectories relaxed to the same side. This result 
indicates that the forward and backward trajectories be- 
have in an almost uncorrelated way as one would expect 
for diffusive barrier crossing [2^ . The average trajectory 
crosses the dividing surface at z = more than 8 times 
and there are trajectories which recross 60 times. This 
large number of recrossings is also indicative of diffusive 
dynamics. 

To characterize the transition mechanism in more de- 
tail, we have analyzed permanence times of the proton 
in the tube after release from z = 0. The distribution of 
these permanence times extracted from 5000 trajectories 
is shown in Fig. HI Here, the permanence time on the bar- 
rier is the time the proton needs to reach |z| = 8 A start- 
ing from the barrier top. The distribution peaks at about 
0.6ps and then decays exponentially with a time constant 
of about G.93ps. This distribution of permanence times 
is reproduced very well by a one-dimensional Brownian 
particle evolving on the effective potential F{z) of Fig. [2] 
with a diffusion constant D = 7 A^ps~^, about half that 
estimated for protons in long water-filled tubes in vac- 
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FIG. 4: Distribution of permanence times of the proton in 
the pore after release from the pore center from simulations 
(circles) and Brownian diffusion on F{z) (solid lines). (Inset; 
logarithmic scale for P{t).) 



uum [8]. The agreement between the molecular dynam- 
ics results for the full system and the one-dimensional 
Brownian dynamics simulation again indicates that the 
proton motion is diffusive and that the center of charge 
is an appropriate reaction coordinate. 

The distribution of permanence times of the proton on 
the free energy barrier can be roughly modeled by a one- 
dimensional diffusion process on a flat potential, starting 
from z = and terminated at ±L/2. The width of the 
almost flat barrier is L « 8 A (Fig. [2]). At long times, 
the resulting distribution of permanence times decays as 
P{t) « ATrDex.p{-TT'^Dt/L'^)/L'^. This exponential decay 
accurately reproduces the long time tail of the distribu- 
tion of permanence times observed in our simulations and 
plotted in Fig. g] 

The resulting picture of a protonic defect diffusing 
through the pore under the influence of the effective po- 
tential F{z) has implications on the design of conducting 
pores. Increasing their leiigth L will reduce the trans- 
mission coefficient as 1/L [27| and hence lower the con- 
ductance. This effect will be enhanced by the larger des- 
olvation penalty arising in longer pores. However, ii sing 
nanopores of higher polarity, possibly embedded 2^, 23| 
in a high-dielectric medium, should greatly reduce the 
desolvation cost and may result in the ideal combination 
of high proton mobility and concentration yielding pro- 
ton current densities comparable to those measured for 
polymer electrolyte membranes. 
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